A simple example of fitting a nonlinear regression and nonlinear mixed-effects model

First we load the NLreg package and and the DataFrames package. The readtable function in creates a DataFrame from a comma-separated or tab-separated file's contents. The file can be compressed, as shown here.



In [11]:

    
using DataFrames,NLreg
const sd1 = readtable(Pkg.dir("NLreg","data","sd1.csv.gz"));









    



Warning: redefining constant sd1

Nonlinear regression models are represented as types in the NLreg package. The type can implement several constructors. Here we generate a BolusSD1 model from a formula expression and a data frame.



In [12]:

    
m = BolusSD1(CONC ~ TIME, sd1);
typeof(m)









    Out[12]:





BolusSD1{Float64} (constructor with 1 method)

This model is a concrete type inheriting from the abstract type PLregModF, a partially-linear regression model function.



In [13]:

    
typeof(m) <: PLregModF









    Out[13]:





true

Add information here about the size of a model, once those methods are defined
Also add formula methods

The formula indicates that the response is CONC and the first (and only) covariate is TIME. The names of the parameters for the model are provided by pnames.



In [14]:

    
pnames(m)









    Out[14]:





2-element Array{ASCIIString,1}:
 "V"
 "K"

The first parameter, V, representing the volume of distribution, is conditionally linear. The second parameter, K, the rate constant is nonlinear.

Fitting a nonlinear regression model

A partially linear least squares fit is represented by the type PLinearLS, which is constructed from a nonlinear regression model and, optionally, initial values for the parameters.

If initial values for the parameters are not given the model's initpars method is used to generate initial values.



In [15]:

    
show(initpars(m))









    



[0.2793108696162248]

Note that the initial estimates are for the nonlinear parameters only.

The fit function fits the model using the Golub-Pereyra algorithm. An optional second argument determines if verbose output is generated.



In [16]:

    
pl = fit(m,true)









    



Iteration:  1, rss = 111.039, cvg = 0.00367759 at [1.20139,0.279311]
   Incr: [-0.0329797]  f = 1.0, rss = 110.597
Iteration:  2, rss = 110.597, cvg = 1.53792e-6 at [1.14412,0.246331]
   Incr: [-0.000604764]  f = 1.0, rss = 110.597
Iteration:  3, rss = 110.597, cvg = 6.08026e-9 at [1.14303,0.245726]
   Incr: [-3.79507e-5]  f = 1.0, rss = 110.597







    Out[16]:





Nonlinear least squares fit to 580 observations

     Estimate Std.Error t value Pr(>|t|)
V     1.14296 0.0495656 23.0595  < eps()
K    0.245688 0.0202414 12.1379  < eps()

Residual sum of squares at estimates: 110.597
Residual standard error = 0.43743 on 578 degrees of freedom

There are several extractor methods such as coef, coeftable, deviance, pnames, stderr and vcov for this type.



In [17]:

    
coeftable(pl) # the coefficient table shown above









    Out[17]:





     Estimate Std.Error t value Pr(>|t|)
V     1.14296 0.0495656 23.0595  < eps()
K    0.245688 0.0202414 12.1379  < eps()



In [18]:

    
deviance(pl)  # residual sum-of-squares









    Out[18]:





110.5972863376179

In this model we used the elimination rate constant, K, as a parameter. Suppose that instead we wish to use the logarithm of the rate constant as a parameter. We could create a new model type or we could combine the existing model with a parameter transformation.

Consider how to code up parameter transformations of both parameters.



In [19]:

    
pl2 = fit([LogTr]*BolusSD1(CONC ~ TIME, sd1))









    Out[19]:





Nonlinear least squares fit to 580 observations

        Estimate Std.Error  t value Pr(>|t|)
V        1.14296 0.0495654  23.0595  < eps()
log(K)   -1.4037 0.0823867 -17.0379  < eps()

Residual sum of squares at estimates: 110.597
Residual standard error = 0.43743 on 578 degrees of freedom

The two fits are equivalent in that they produce the same sum of squared residuals and the parameter estimates can be transformed from one scale to the other.



In [20]:

    
using Base.Test
@test_approx_eq_eps coef(pl) [coef(pl2)[1], exp(coef(pl2)[2])] 1e-5

Fitting a simple nonlinear mixed-effects model

A simple nonlinear mixed-effects model has a random effect for each parameter for each group (e.g. subject). Uncorrelated random effects are indicated by specifying a Diagonal $\Lambda$ matrix. Correlated random effects user a Triangular $\Lambda$ matrix. The model form and initial estimates of the fixed-effects parameters can be taken from a fitted NonlinearLS object.

use an extractor function for mmf



In [21]:

    
nm = fit(SimpleNLMM(pl2,array(sd1[:ID]),Diagonal))









    



type CompositePLModF has no field mmf
while loading In[21], in expression starting on line 1
 in updtmu! at /home/bates/.julia/NLreg/src/nonlinreg.jl:103
 in prss! at /home/bates/.julia/NLreg/src/simplenlmm.jl:54
 in pnls! at /home/bates/.julia/NLreg/src/simplenlmm.jl:100
 in pnls! at /home/bates/.julia/NLreg/src/simplenlmm.jl:117
 in fit at /home/bates/.julia/NLreg/src/nlmm.jl:51



In [22]:

    
nm1 = fit(SimpleNLMM(pl2,array(sd1[:ID]),Triangular))









    



type Triangular is immutable
while loading In[22], in expression starting on line 1
 in SimpleNLMM at /home/bates/.julia/NLreg/src/simplenlmm.jl:37

Although not currently shown in the output, the within-subject correlation of the random effects for this model fit is -0.091. A likelihood ratio test comparting nm versus nm1 has a p-value of 39.6% from which we would conclude that the more complex model, nm1, does not provide a sufficiently better fit to justify the additional parameter.